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Abstract 

In this paper, we present a novel Lagrangian formulation of the equations 
of motion for point vortices on the unit 2-sphere. We show first that no linear 
Lagrangian formulation exists directly on the 2-sphere but that a Lagrangian 
may be constructed by pulling back the dynamics to the 3-sphere by means 
of the Hopf fibration. We then use the isomorphism of the 3-sphere with the 
Lie group SU (2) to derive a variational Lie group integrator for point vortices 
which is symplectic, second-order, and preserves the unit-length constraint. At 
the end of the paper, we compare our integrator with classical fourth-order 
Runge-Kutta, the second-order midpoint method, and a standard Lie group 
Munthe-Kaas method. 

1 Introduction 

Point vortices are point-like singularities in the vorticity field of an ideal fluid. First 
described by von Helmholtz [1858], they form a finite-dimensional singular solution 
of the Euler equations and are now a classical subject in hydrodynamics, see among 
others Lamb [1945]; Milne-Thomson [1968]; Saffman [1992]; Newton [2001]. 

The interest in point vortices is two- fold. On the one hand, paraphrasing Aref [2007], 
the description of point vortices forms a veritable playground for classical mathe- 
matics and gives rise to interesting phenomena from dynamical systems, such as 
periodic motions (Souliere and Tokieda [2002]; Borisov, Mamaev, and Kilin [2004]), 
(relative) equilibria (Kidambi and Newton [1998]; Aref [2011]), and chaotic advec- 
tion and topological chaos in fluids (Boyland, Stremler, and Aref [2003]). On the 
numerical front, on the other hand, desingularizations of the point vortex equa- 
tions, such as the classical vortex blob method of Chorin [1973] form the basis for 
important classes of particle methods for the Euler and Navier-Stokes equations. 
The idea is that the vorticity field of an arbitrary fluid can be approximated by a 
number of vortex blobs whose motion is then followed in time. Strong analytical 
estimates exist that link the behavior of the vortex blobs with the solution of the 
Euler equations that they approximate (Majda and Bertozzi [2002]). 
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On the sphere, the dynamics of point vortices was first described by Bogomolov 
[1977] after a model by Gromeka (see Newton [2001] for an historical overview) and 
is in some sense a generalization of the planar case (see also Kimura and Okamoto 
[1987] and Polvani and Dritschel [1993]). The relevance of point vortices of the 
sphere lies in the fact that they provide a first approximation of the behavior of 
certain geophysical flows for which the curvature of the earth is important, and 
which persist over long periods of time. The mathematical description of point 
vortices on the sphere is an area of active research: fixed and relative equilibria of 
the three-vortex problem were described in Kidambi and Newton [1998] (see also 
Pekarsky and Marsden [1998]), while more general equilibria were described in Lim, 
Montaldi, and Roberts [2001]; Chamoun, Kanso, and Newton [2009]; Newton and 
Sakajo [2011]. Conditions for the collapse of point vortex configurations on the 
sphere were established in Kidambi and Newton [1998] and Sakajo [2008]. 

Most of the research on point vortices on the sphere has focused on the existence 
of analytical solutions such as relative equilibria for few point vortices, but com- 
paratively little is known about the behavior of arbitrary configurations of vortices. 
One of the contributions of this paper is to construct a geometric numerical inte- 
grator which is second-order accurate, preserves the geometry of the sphere, and is 
symplectic. As symplectic integrators are known to capture the long-term behavior 
of a Hamiltonian system better than classical integrators (see Marsden and West 
[2001]; Hairer, Lubich, and Wanner [2002]), we expect our geometric integrator to 
give insight into the behavior of non-equilibrium vortex configurations, even over 
long integration times. 

1.1 Aims and contributions of this paper 

The contributions of this paper are two-fold. In the first part of this paper, we 
construct a Lagrangian description for point vortices on the sphere in terms of pairs 
of complex numbers. We first review the Lagrangian description for point vortices 
in the plane (see e.g. Chapman [1978]; Newton [2001]; Rowley and Marsden [2002]) 
and then show via a simple topological argument that no (linear) Lagrangian exists 
for the dynamics of point vortices on the two-dimensional sphere S 2 . 

We then use the Hopf fibration, a distinguished submersion from the three-sphere 
§ 3 to the two-sphere § 2 , to pull back the Hamiltonian description to § 3 , where the 
topological obstruction for the existence of a linear Lagrangian vanishes. We explic- 
itly construct this Lagrangian and we show that the equations of motion give rise to 
a (finite-dimensional) non-linear Schrddinger equation on § 3 with gauge freedom. 
These equations bear a remarkable similarity to the equations of motion for point 
vortices in the complex plane, only now the location of each point vortex is specified 
by a pair of complex numbers (or equivalently, a (unit) quaternion) instead of a 
single one. 

In the second part of the paper, we design a variational numerical integrator for point 
vortices on the sphere using the linear Lagrangian on S 3 . We use the identification 
between the 3-sphere § 3 and the Lie group SU(2) of special unitary 2-by-2 matrices 
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to write the update equation for the integrator as a fixed-point equation in the Lie 
algebra su(2), and we show how the discrete equations of motion are symplectic, 
self-adjoint, second-order, and preserve the unit-length constraint in S 3 . At the 
end of the paper, we compare our integrator to the classical 4th order Runge-Kutta 
method, as well as to a number of geometric integration methods. We show that the 
geometric integrators, and in particular the Hopf variational integrator, outperform 
Runge-Kutta in the medium run, even though they are only second-order accurate. 



1.2 Background and historical overview 

Linear Lagrangian formulation for planar vortices. We review here the La- 
grangian and Hamiltonian descriptions of point vortices in the plane. The Hamilto- 
nian description of point vortices on the sphere will be reviewed in Section 3. 

For point vortices in the plane, the equations of motion are given in complex form 

by 

Here, the z a (a = 1, . . . , N) represent the locations in the complex plane of each of 
the vortices, and T a is a real parameter which specifies the circulation around each 
vortex. The Hamiltonian function is given by 

H(zi, ...,z N ) = -— y~] T a T/3 log \z a - zp\ 2 . (1.2) 

a</3 

These equations can be derived from a Lagrangian which is linear in the velocities 
(see Chapman [1978]) and is given by 



1 N 

L = — ^2 ^ a {z* a z a - z a z* a ) - H(zi, ...,z N ). (1.3) 

a=l 

For future reference, we point out that the linear part of the Lagrangian can be 
written as ^ T a 6(z a , z a ), where 9 is the one-form given by 

9 = ^lm(z*dz). (1.4) 

The exterior derivative of 6 is nothing but the area form on the complex plane: 

d8 = - lm(dz* A dz) = dx A dy. (1.5) 

It can be shown that the flow of the point vortex equations (1.1) preserves a weighted 
sum of such area forms, given by 

N 

} y F a d9 a = T a dx a A dy a , 

a=l a=l 
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where d9 a refers to the area form (1.5) expressed in the coordinates of the ath 
vortex. This is an example of a symplectic form on the phase space C . 

The advantage of having a Lagrangian description for the dynamics of point vortices 
is that the standard results for the construction of Lagrangian variational integrators 
(see Marsden and West [2001] for an overview) can now be applied. This is the key 
observation of Rowley and Marsden [2002] , who constructed a class of second-order 
variational integrators by discretizing the Lagrangian (1.3) using centered finite 
differences. 

Before turning to the case of point vortices on the sphere, we point out that many 
non-canonical Hamiltonian systems can be rephrased as Euler-Lagrange equations 
that come from a Lagrangian which is linear in the velocities. This observation was 
made by Birkhoff [1966] in his study of Pfaffian systems and was used in Faddeev and 
Jackiw [1988] as a starting point for the description of Hamiltonian systems with 
constraints. Linear Lagrangians also appear in the description of the non-linear 
Schrbdinger equation and the KdV equation. 



The dynamics of point vortices on the sphere. The equations of motion for 
iV point vortices with strengths Ti, i = 1 , . . . , N on the unit sphere § 2 can be written 
as follows (see Newton [2001]). If we denote the position vector of the ith vortex by 
Xj (so that ||xj|| = 1), the point vortex equations can be written in Euclidian form 
as 

x * = ^X> 1 + cj2 _ Xfe . Xj , (1-6) 

where a is a small regularization parameter which is added to ensure that the limit 
of the right-hand side exists when x^ tends to Xj. Note that the equations (1.6) 
conserve the vortex moment, defined as 

N 



i=l 



Non-existence of a linear Lagrangian for vortices on the sphere. In Sec- 
tion 3, we will review the Hamiltonian formulation for the point vortex equations 
(1.6). We now discuss the Lagrangian formulation, and in particular we argue that 
no linear Lagrangian exists for the dynamics of point vortices on § 2 . 

For point vortices on the sphere, the phase space is the product (§ 2 ) N of N copies 
of the unit sphere § 2 , equipped with a symplectic form which is a weighted sum 
of the area forms on the individual spheres. This particular symplectic form is not 
exact, due to S 2 being compact; the standard way to see this is that if it were, then 
O, = d© for a one-form on the sphere, and therefore Stokes' theorem would give 

4vr= / n= [ de= / 9 = 0, 



a contradiction. 
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This rules out the existence of a linear Lagrangian on 8 2 , as such a Lagrangian 
would necessarily have to be of the form L = — H, with O so that dO is the 
symplectic form. However, we will see that the symplectic form Q on § 2 can be 
pulled back to the three-sphere § 3 so that the resulting form is exact. This will 
allow us to construct a linear Lagrangian for point vortices on § 3 . By discretizing the 
Lagrangian variational principle on § 3 (see Marsden and West [2001]), we will then 
be able to construct a variational integrator for point vortices which is automatically 
symplectic, second-order, and unit-length preserving. 

Other approaches to the numerical integration of point vortices. Hamilto- 
nian variational principles have been developed by Oh [1997] in the context of Floer 
homology and by Novikov [1982] for Morse theory (see also Cendra and Marsden 
[1987]). On the numerical front, geometrical numerical integration of Hamiltonian 
systems was described in Brown [2006], Ma and Rowley [2010] and Leok and Zhang 
[2011], but all of these references assume that the underlying symplectic manifold is 
exact. For non-exact symplectic forms (e.g. the case of point vortices on the sphere) 
it is as of yet not clear how to discretize the Hamiltonian variational principle so that 
the resulting numerical algorithms share some of the properties of the continuous 
system (such as symplecticity and momentum preservation). We do remark that 
Ma and Rowley [2010] perform a similar pullback as in this paper, but using the Lie 
algebra of the rotation group SO (3) instead of the special unitary group SU(2), in 
order to make the dynamics of point vortices on the sphere amenable to geometric 
integration. 

2 The Hopf fibration 

In this section, we introduce our notation and review some aspects of the geometry 
of the spheres § 2 ,§ 3 and the Hopf fibration. This material is standard and can 
be found in any geometric physics textbook, for instance Frankel [2004]. More 
information about the Hopf fibration and its role in physics and geometry can be 
found in Montgomery [2002]; Urbantke [2003]; Lemaitre [1948] and the references 
therein. 

Notation. We will denote vectors in C 2 and their Hermitian conjugates by 

and (p' := [z* , it*] , 

where z* is the complex conjugate of z G C. The Hermitian conjugate of a complex 
matrix A will be denoted by A'. 

Lowercase Roman letters a, b, . . . will refer to the components <p a of a vector ip in 
C 2 . The Greek letters a, (3, . . . will refer to the Cartesian components x a of a vector 
x G R 3 . The imaginary unit will be denoted by i. 
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The Hermitian inner product on C 2 is given by 

(<Pi, <P2) ■= ^1^2 = z\z 2 + u\u 2 . 
Note that the Euclidian inner product on C 2 can be expressed as 

Re (ipi, <p 2 ) = Rb{z{z 2 +u*u 2 ). (2.1) 

The geometry of S 2 . We think of the two-sphere § 2 as the set of all unit-length 
vectors x in M 3 . The tangent plane T x § 2 at an element x 6 S 2 consists of all vectors, 
denoted by ix G M 3 , which are orthogonal to x: 

T x § 2 = {5x £ M 3 : x- 5x = 0}. 

In Cartesian coordinates, the area form f2 on S 2 can be described as follows: Q is 
the differential two-form given by 

fi(x)(£x,<5y)=x-(<5xx<5y), (2.2) 

for all x£§ 2 and <5x, 5y G T x § 2 . In spherical coordinates, O = s'mOdO A dep. Note 
that $7 is not exact. 



The geometry of § 3 and the group SU{2). We let S 3 be the unit sphere in 



S 3 = {(z, u) g 



I I 2 i I |2 

\z\ + U 



!}■ 



The tangent plane at an element (p £ S 3 is given by the set of all vectors, denoted 
by 5f £ C 2 , which are orthogonal to cp: 



{Sip £ C 2 : Re{S(p,ip) = 0}, 



(2.3) 



where we have expressed the orthogonality between <p and dip using the inner product 
(2.1) in C 2 . 

The unit sphere S 3 can be embedded into the complex 2-by-2 matrices, by means of 
the mapping 



-u, 



€ M 2 (C), 



whose range is precisely the Lie group SU(2) consisting of all Hermitian matrices 
([/• = U) with unit determinant (det U = 1). The Lie algebra of SU(2) is the vector 
space su(2), consisting of all 2-by-2 matrices A which are anti- Hermitian {A^ = — A) 
and traceless (trA = 0). The identification of S 3 with SU(2) provides a convenient 
description for the tangent spaces (2.3): we have that 5ip £ T^S 3 if and only if there 
is a matrix A £ su(2) such that 

5ip = Aip. (2.4) 

To see this, note that A^ = —A implies that (ip, Aip) is purely imaginary, so that 
Re ((p, Aip) = 0. 



2 The Hopf fibration 7 



The Lie algebra su(2) has dimension 3 and a convenient basis is given by the matrices 
= 1,2,3, where the a Q are the Pauli spin matrices: 





"0 


f 




"0 -i" 


(71 = 


1 





, 0-2 = 


i 



and c 3 



1 

-1 



(2.5) 



Given a matrix A £ su(2), we will denote its components in this basis by a a , 
a = 1, 2, 3, and we put a := (ai, 02, 03). Explicitly, 



A = a • (icr) = a a l 



10"n 



(2.6) 



0=1 



where a represents the vector (a 1 , a 2 , a 3 ). We will refer to a 6 M 3 as the vector 
representation of the matrix A £ su(2). 

The Pauli matrices satisfy a number of useful identities: the multiplication identity 



is 



a a ap = 5 a pl + i e a ^a. 

7=1 



71 



(2.7) 



for a, (3 = 1, 2, 3, where I is the 2-by-2 unit matrix and e Qj g 7 the Levi-Civita permu- 
tation symbol. Secondly, there is the completeness property 



"^^{o~a)ab{o~a)cd = ^adhc ~ Sabred, 



(2.8 



a=l 



for all a,b,c,d = 1,2. Proofs of these identities can be found in any standard 
textbook on quantum mechanics. 



The Hopf fibration. The group U{1) = S 1 of unit complex numbers acts on S 3 
by the diagonal action: e %e ■ (z,u) = (e l6 z, e l6 u) for all e ld S S 1 and (z,u) G S 3 . In 
terms of SU (2)-matrices, this action is described as 



z —u* 




z 


-u*~ 




~e w " 


* 

u z 


■e i8 = 


u 


z* 




e~ w 



(2.9) 



The orbit space of this action, S 3 /S 1 , can be identified with the two-sphere § 2 . 
Explicitly, there exists a surjective submersion it : S 3 — > § 2 , called the Hopf fibration, 
given by 

tt(z,u) = (2Re(z*u),2lm(z*u),\z\ 2 - \u\ 2 ), (2.10) 

and the fibers of tt coincide with the orbits of the group S 1 in §> 3 . In geometrical 
terms, the map ir : S 3 — > E 2 makes S 3 into the total space of a right principal fiber 
bundle with structure group S 1 over §> 2 . We will refer to the orbits of the S 1 -action 
(2.9) as the S 1 -fibers of S 3 . 
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The Hopf map can be expressed conveniently in terms of the Pauli matrices as 
follows. We let a be the vector (o"i, 02, 03). The Hopf map (2.10) can then be 
described as 

ir(tp) = <p*o-<p. (2.11) 

The right-hand side should be interpreted as a vector x in i 3 , whose components 
are given by x a := (p'a a (p, a = 1, 2, 3. 

The inner product of two vectors x, y £ M 3 can be given a convenient description 
using the Hopf map. Let x = ip^trip and y = tp^crip. A straightforward consequence 
of (2.8) is then that 

x • y = 2(pty)tyV) - (pVXV'ty)- (2-12) 

Connection one- form and curvature. On § 3 , there is a distinguished one- form 
A which will play a crucial role in obtaining the Lagrangian formulation for point 
vortices. Explicitly, it is given by 

A((p) = Im(^ dip), 

and we denote the contraction of A(ip) with a vector <jp = (i, u) by 

A{<p) ■ <p = Im((^V) = Im(z*i + u*u). (2.13) 

Note the similarity between A and the one- form 6 given in (1.4). 

The form A is the connection one-form of a principal fiber bundle connection, but 
we will just treat it as a one-form. The curvature of A is given by 

d„4 = Im(oV t A dip) = lm(dz* A dz + du* A du), 

and it can be shown that the area form Q on § 2 , given by (2.2), satisfies 

ir*Q = 2dA. (2.14) 

This result states that the two- form f2, which is not exact, nevertheless becomes 
exact when pulled back along the Hopf map to S 3 . This will allow us to construct 
a linear Lagrangian for point vortices on S 3 . 



3 Hamiltonian formulation of the vortex equations 

In this section, we review the Hamiltonian description of the equations of motion 
(1.6) for point vortices on the unit sphere. This system of first-order ODEs can be 
written in Hamiltonian form, where the phase space is the Cartesian product (S 2 ) N 
of N copies of the unit sphere § 2 , equipped with the symplectic form 



N 

#r(xi, . . . ,Xjv) = y^T^(x^), 
i=i 



(3.1) 
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where f2 is the standard symplectic area form on § 2 , given by (2.2). 
The Hamiltonian function is given by 

^ = -^E r ^ lo g( 2cj2 + ^)' ( 3 - 2 ) 

i<j 

where Uj := || the chord distance between the ith and the jth vortex and 

a is the cutoff parameter introduced in (1.6). 

Hamilton's equations are then given by ix^r = dH. Explicitly, we are looking for 
a curve t h-» x(t) G (§ 2 ) N such that, for any variation <5x(i) € T x ( t \(S 2 ) , we have 
that 

£ r (x,5x) = dH(x)-<5x. 
Using the expression (2.2) for the symplectic form, this can be written as 

N N 
i=l i=l 

so that 

I\(xi x ii) = V XiJ ff(x) + A iXi , (3.3) 

where the Lagrange multipliers Aj, i = 1, . . . , N, have been introduced to enforce the 
constraint that the variations <5xj be tangent to the unit sphere, so that Xj • <5xj = 
for alH = 1, . . . , N. Taking the cross product of (3.3) with x» then results in 

Vi(xi x ±i) x x, = V Xi H(x) x Xi, 

and after expanding the double cross product and using the fact that ||x$|| = 1, we 
obtain 

Ti±i = V XiJ ff(x) x Xi , (3.4) 

which is equivalent to (1.6). 

4 Lagrangian formulation of the vortex equations on § 3 

In this section, we show how the Hamiltonian equations (1.6) for point vortices can 
be given a Lagrangian formulation. To do this, we lift the point vortex system from 
the two-sphere S 2 to the three-sphere § 3 using the Hopf fibration. 

Pullback of the Hamiltonian H. Using the projection ir given in (2.10), we 
may pull back the Hamiltonian function on S 2 to S 3 . If we denote the Hamiltonian 
function (3.2) by H§2 and the pullback by H§3, then we have that H§a = ir*H§2, or 
explicitly, 

H s3 (tpi,...,ip N ) = H S 2(w((pi), . . . ,it(ipn)), (4.1) 
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for all (pi, . . . , LfN € § 3 . A straightforward computation shows that if §3 is given by 

flss (pi,...,^) :=-^XJ r ^ lo s[ 2<72 + 4 ( 1 -|^'^')| 2 )] • ( 4 ' 2 ) 

In the remainder, we will drop the subscript '8 3 ' on the Hamiltonian function, 
denoting H§3 simply as H. Note that H is invariant under multiplication by e l9 G S 1 
in each argument separately: 

H(...,j <p k ,...)=H(...,<p k ,...), (4.3) 

for k = 1, . . . , N. The infinitesimal version of this symmetry is 

a^ n -' fil M = ' (4 ' 4) 

where there is no sum over the index k. 

Since multiplying ip^ by a phase factor e 10 corresponds to moving along the S^fiber 
through ipk, we have that H depends only on the chord distance between the S 1 - 
fibers through <pi, . . . , tpjsf. This can be shown explicitly in (4.2) by expressing the 
inner product \(ifi, (pj)\ in terms of the Euclidian distance D((pi,(fj) between the 
S 1 -fibers through ipi and ipj, where 

T>(fi,fj) = 2(1 - \((pu(pj)\). 



The linear Lagrangian and the equations of motion. We now have all the 
elements to formulate a Lagrangian description for point vortices on S 3 . Recall that 
a linear Lagrangian has the general form L = O — H, where dO is the symplectic 
form. The symplectic structure on (S 3 )^ is given by the pullback of the symplectic 
structure on (S 2 )^, 



so it follows that = 2 5^i=i ^iAi. Therefore, we obtain 

N 

L = 2 FiA{ipi) ■ (pi - H(<pi, <p N ), (4.5) 
i=i 

where ipi G § 3 for i = 1, . . . , N. This generalizes the expression (1.3) for the linear 
Lagrangian for point vortices in the plane. 

The action functional is defined as 

«%>(•))= L(<p(t),<p(t))dt, (4.6) 
J t 



4 Lagrangian formulation of the vortex equations on § 3 
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where ip(t) := (v?i(i), . • . ,(pw(t)) is a curve in (S ) defined on the interval [io^iL 
and its variation is given explicitly by 



i= i V ^,7 »=i 



(4.7) 



where the infinitesimal variations 5<£>j and <5(/?J need to be prescribed carefully. Since 
<Pi is an element of § 3 , the variations Sipi are elements of T^S 3 . Specifically, we have 
that S(fi is orthogonal to </?j. This relation may be incorporated using Lagrange 
multipliers Aj, resulting in the Euler-Lagrange equations 



2'iTiipi 



dH 



together with their Hermitian conjugates and the unit-length constraints 



1. 



(4- 



(4.9) 



This equation is the analogue of (1.1) for vortices on § 3 and can be seen as a nonlinear 
Schrddinger equation on the product space (S 3 )^. The analogy with (1.1) can be 
made more striking by interpreting ifi as a unit quaternion, so that (4.8) becomes 
(up to a constant) the quaternionic version of the complex equation (1.1). 



Determining the Lagrange multipliers. A curious feature of these equations 
is that the multipliers Aj reflect gauge degrees of freedom, that is, any choice of Aj 
will preserve the unit length constraint equally well. To see this, take the time 
derivative of (4.9) and substitute the equations of motion; the result is 




which simplifies to 



dH t OH n 

d<Pi dip] 



from which Aj is absent. This expression is nothing but the infinitesimal symmetry 
relation (4.4) and is therefore identically satisfied. 

With hindsight, it is not surprising that there is some indeterminacy in the solutions 
of (4.8). After all, these equations arise as pullbacks of equations on S 2 . From this 
point of view, changing the multipliers Aj will change the dynamics in the direction 
of the S^fibers, but will leave the horizontal dynamics (which projects down to § 2 ) 
unchanged. This is similar to the un-reduction approach of Bruveris, Ellis, Holm, 
and Gay-Balmaz [2011], in which a complicated dynamical system on a manifold X 
is mapped into a simpler problem on the total space of a principal fiber bundle over 
X. 



4 Lagrangian formulation of the vortex equations on § 3 
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Relation with the equations on (E> 2 ) N . By construction, the flow of the equa- 
tions (4.8) on (§> 3 ) N will project down onto the flow of the point vortex equations 
(1.6) on (S 2 ) N . It is instructive, however, to see this explicitly. 

We start again from the variational principle (4.7), but now we do not introduce 
a Lagrange multiplier to incorporate the unit-length constraint. For the sake of 
clarity, we suppress the explicit index i in (4.7) to arrive at 



5S = Sip 



- 2ir ^ + vJ + v 



(4.10) 



As the variation dip is tangent to § 3 , it can be written as Sip = Aip, where A G su(2); 
see (2.4). Similarly, we have that 5^ = ip'A' = —ip^A. Upon substituting these 
expressions in (4.10), we arrive at 

« = Vi»(-ar* + ^) + (!iir*' + ^)^ 



2 Re 



so that 5S = for all A G su(2) if and only if 





7 ^ t 9H\ 
\ 2lTip ^ 




Re 


= 0, a = 1,2, 3, 



(4.11) 



where the a a are the Pauli matrices (2.5). Note that these equations are equivalent 
to (4.8). 

We now let x G S 2 be the image of (p G S 3 under the Hopf map, and we recall 
from (2.11) that the components of x are given by x a = p^a a ip. Taking the time 
derivative, we obtain 



x a = 2 Re {^a a p 



(4.12) 



Similarly, we recall that the Hamiltonian functions Hg2 and H S 3 are related by 
(4.1), or explicitly #§3(92,^) = H§2{<p^a a <p). Taking the derivative with respect to 
ip^ yields 

~W = -tof**"' (413) 

and a small calculation, involving the multiplication identity (2.7), then shows that 



Re 



dH 



\ - dH S 2 

^ taPl dx Xl 

13,7 P 



(V x #§2 X x) a . 



Substituting this expression and (4.12) into (4.11) then results in the following vector 
equations: Tx = V x f/§2 x x, which, upon restoring the sum over all vortices, are 
nothing but the point vortex equations (3.4). 
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Intrinsic formulation of the point vortex equations. The point vortex equa- 
tions (4.8) on (S 3 )^ can be written in geometric form as follows. The point vortex 
trajectories are integral curves of a vector field T on (S 3 ) N , determined implicitly 
from the equation 

i T B$ 3 =dH S3 , (4-14) 

where B§ 3 is the pullback to (S 3 )^ of the point vortex symplectic form B given 
in (3.1): B§ 3 = ir*B, with ir : S 3 — > E 2 the Hopf map. This makes it clear that 
the equations (4.8) do not determine the dynamics completely: we may add to 
the solution T of (4.14) any ir- vertical vector field without changing the projected 
dynamics. 

We may resolve this indeterminacy by replacing S 3 by its symplectification, which 
is the product manifold S 3 x M + equipped with the symplectic form 

B = d(rA), 

where r is the coordinate on the M + -factor. The motion of point vortices on this 
twice-enlarged space projects down to the motion of point vortices on S 2 , and can be 
viewed, paraphrasing the terminology of Kostant [2005], as aversion of "prequantum 
vortex dynamics." 1 



5 Variational integrators on SU(2) 

In this section, we propose a discrete version of the point vortex equations on (S 3 )^. 
We begin by discretizing the linear Lagrangian (4.5) using centered finite differences. 
By taking discrete variations, we then obtain a discrete version of the equations (4.8) 
where the constraints are enforced using a Lagrange multiplier. These equations can 
be seen as a version of the Moser-Veselov equations (see Moser and Veselov [1991]) 
on (S 3 )^. 

By projecting onto the annihilator space of the constraint forces, we then obtain a 
discrete version of the projected equations (4.11). Finally, we use the isomorphism 
between § 3 and SU (2) to write the discrete equations of motion in the form of a 
Lie group variational integrator (see Lee, Leok, and McClamroch [2007]) on SU(2) 
and we argue that this form of the equations is especially well-suited for numerical 
implementation. 

5.1 Discrete Lagrangian and discrete equations of motion 

Let M be the number of discrete time steps, with constant time increment h > 0, 
and denote the variables at time t n := nh by (p n := (99™, . . . , ip 1 ^) S (S 3 )^. We now 
propose a discrete counterpart of the linear Lagrangian L in (4.5) by approximat- 
ing the action functional (4.6) over the interval [t ra ,i ra +i] by using piecewise linear 



1 We thank M. Gotay for bringing this point to our attention. 
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interpolants and the midpoint rule (see Marsden and West [2001]). In this way, we 
construct a discrete Lagrangian L d : (S 3 )^ x (S 3 )^ — > K of the form 

/ml J.,/,1+1 in n+l — tn n \ 

W, v n+1 ) = hL , * h V ) ■ (5.1) 

Explicitly, the discrete linear Lagrangian is given by 

JV 

L d (<p n ,<p»+l) =2Y j r i A(tf +1/2 ) • « +1 - + hH{^ n+1 / 2 ), (5.2) 
i=l 

where (^ n+1 / 2 := (<^ n + <p n+1 )/2. 

The discrete action sum is now defined as 

M-l 

i=i 

and taking variations with respect to ipf and (iffy yields 
N M-l 

*s* = £ E [ - ir * - ^r 1 ) - hxn ^ 

k=l n=l 

+ £ (d^H{^-^) + DjHitf+W)) } + (ex.), (5.3) 

where "(c.c.)" stands for the complex conjugate of the previous expression. The 
discrete Euler-Lagrange equations are hence 



-IT* (^ +1 - ^r 1 ) + \ (D^Hi^- 1 ' 2 ) + D A H{v n+l/2 )) ~ fcAM = 0, 



(5.4) 



together with their Hermitian conjugates and the unit-length constraints 

(^ +1 ,< +1 > = 1, (5.5) 

and can be viewed as the discrete analogues of the continuous equations (4.8). 

In contrast to the continuous case, the Lagrange multipliers in (5.4) are no 
longer arbitrary. Instead, they can be found by substituting the discrete equations 
of motion into the unit-length constraint (5.5) and solving the resulting system of 
quadratic equations for A£. 



Other discretizations. Instead of the midpoint quadrature formula, leading to 
the discrete Lagrangian (5.1), we could have chosen another approximation for the 
discrete Lagrangian, leading to a different set of discrete equations. Although these 
equations formally exhibit the same properties (symplecticity, preservation of the 
unit length constraint, etc.) as the midpoint method introduced previously, some 
of them suffer from undesirable side-effects. A particularly revealing example is 
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obtained by using the trapezoid rule instead of the midpoint rule, leading to the 
discrete Lagrangian 



2 V V h J V h 

which gives rise to the discrete Euler-Lagrange equations 

- ir fc - pT 1 ) + hD^H(^) = hXttf, (5.6) 

together with unit-length constraints (5.5). This system, however, is equivalent to 
the 2-step symmetric explicit midpoint integrator (see Hairer, Lubich, and Wanner 
[2002, Sec. XV. 3. 2]), which is well-known to exhibit parasitic solutions which grow 
linearly in time. The same observation is true for the integrator of Rowley and 
Marsden [2002] with a = or 1; see appendix A for details. 



5.2 The projected midpoint equations 

It is possible to obtain an equivalent version of the equations (5.4) which does 
not involve the Lagrange multipliers To this end, we return to the discrete 
variational principle (5.3). As in the case of the continuous variational principle, we 
again suppress the vortex index k to arrive at 

M-i , 

5s d = ^ %»)t [-ir (^ +1 - ^)+ \ (iv^%"~ 1/2 ) + (¥> n+1/2 )) ] +( c - c - 

n=l 

Instead of introducing a Lagrange multiplier to enforce the unit- length constraint, 
we impose the condition that the variations 5ip and Sip' are tangent to the sphere, 
so that they can be written as 

5ip = Aip, and 8<p' = —ip^A, 

where A E su(2) is arbitrary; see (2.4). Proceeding as in the case of the continuous 
variational principle, we then arrive at the following discrete equations: 



Re 



(<P n )Hi* a ) (-ir {Lp n+1 - v n ~ l ) + \ (d^h^-v 2 ) + D^H(<p n +y 2 ))\ 



(5.7) 

for a = 1,2, 3. These equations are supplemented by the unit-length constraint (5.5). 
Note that the equations (5.7) are the discrete version of the continuous projected 
vortex equations (4.11). Another, more direct way to arrive at these equations is 
simply to project the discrete Euler-Lagrange equations (5.4) onto the subspace 
orthogonal to (p n . 



2 Here a refers to the interpolation parameter used in Rowley and Marsden [2002], and should 
not be confused with the cutoff parameter used in the rest of the current paper. 
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The equivalent first-order system. The equations (5.7) are second-order dis- 
crete equations: given (ip 11 ^ 1 , tp n ), the equations can be solved to find ip n+1 . The 
discrete equations of motion are hence not self-starting: given the initial positions 
ip° for the point vortices, a standard one-step integrator needs to be used to find the 
positions ip 1 at the intermediate time t\ = h. Afterwards, the discrete equations of 
motion can be used to integrate the system forwards in time. 

We now present a way to circumvent this problem by writing the two-step method 
(5.7) as the composition of two one-step methods, which are obviously self-starting. 
That this decomposition is possible is a consequence of the fact that the Lagrangian 
(4.5) is linear in the velocities, and will be analyzed further in Appendix A. For 
now, we just focus on the computations for the point vortex equations. 

We first write the equations (5.7) as 



Re 



(V' 



-ilW 



n+1 



n+1/2* 



(5.8) 



where depends on the configurations (f n and ip n only: 



dl := Re 



(5.9) 



We now focus on the left-hand side of (5.8), which we solve for ip n+l : we introduce 
a map <&h '■ § 3 § 3 defined by the property that ^f l (ip n ) = ip n+1 if and only if the 
left-hand side of (5.8) vanishes: 



Re 



(^)t(ia Q ) (-ir(^ +1 - <p n ) + ^ t ^% n+1/2 )) 



0. 



(5.10) 



Now, the adjoint of <3?/i is the map <5>* h : S 3 -»• S 3 defined by 



if and only if <&^f l { i p r ' 



_n-l 



A small calculation then shows that &* h is implicitly determined by the left-hand 
side of (5.9). In other words, 5>t(y> n ~ ) = <p n if and only if 



Re 



0. 



(5.11) 



A possible way of solving the second-order equations (5.7) is therefore as follows. 
Starting from an arbitrary element tp^ 1 £ S 3 we apply alternatingly the adjoint 
map to obtain (p n = ^* h ((p n ~ 1 ), and the direct map $>h to obtain ip n+1 = $/j((/? n ). 
By the previous calculation, the resulting triple (ip n ~ 1 , ip n , ip n+1 ) is then a solution 
of (5.7). We summarize this in the following proposition. 

Proposition 5.1. Let ip 11 ^ 1 G S 3 and put ip n = $^((/9 n_1 ) and ip n+1 = ® h ((p n ). 
The triple (ip n ~ 1 , ip n , <p n+1 ) is then a solution of the discrete second-order equations 
(5.7). 



5.3 Implementing the unit-length constraint: the Cayley map 17 



The maps &h and both approximate the continuous flow, but it can be checked 
that they are only first-order accurate. It is only their composition 3^ o $? that is 
symmetric and hence second-order accurate. 

Note also that this decomposition gives us only a subset of all the solutions of the 
discrete second-order equations (5.7). For instance, if we start from (ip n ~ 1 ,Lp n ) £ 
S 3 x S 3 such that ip n is only approximately equal to Q* h {<p n ~ l ), then (^ n_1 , <f n ) will 
satisfy (5.9) for a non-zero value of d™. However, we can then solve (5.8) for ip n+1 
using that (f™ as the right-hand side. The constants <i™ can now be incorporated in 
the definition of the first-order maps ^h,^, so that solving (5.7) is equivalent to 
taking one step with the first-order map and one step with its adjoint. From this 
point of view, d™ becomes a slack variable. In practice, this will not be an issue, 
since we mostly just start from the initial positions (p° of the point vortices and then 
solve for ip l using (5.11), so that d™ = 0. 



5.3 Implementing the unit-length constraint: the Cayley map 

When we solve (5.10) for f n+1 , we still need to impose the unit-length constraint 
(5.5). This can be done conveniently using the geometry of 577(2): we write the 
update map ip n \-t ip n+1 as 

ip n+1 = U {n) ip n , 

where is an element of SU (2). This ensures that the length of (p n stays constant 
over time, since 

((Z^ 1 ) V+ 1 = ( (/3 n )t([/(«))t[/W^™ = (<p n )i(p n , 

so that, in particular, \\ip n \\ = 1 implies that ||(/2 n+1 || = 1. 
The equations (5.10) for ip n+1 can now be expressed as 



Re 



(^)t(kx Q ) (-T(UW - I 2x2 )<P n + ^IV/^ n+1/2 )) 



0. (5.12) 



These equations can be solved for directly, but a computationally more advan- 
tageous approach is as follows. As long as the step size h is small, the update matrix 
will be in a neighborhood of the identity element in SU(2). We now parametrize 
that neighborhood by means of the Cayley transform Cay : su(2) — > SU(2), given 
by 



Cay(A) = (I + A)(I - A) 



-i 



That is, we search for an element A^ G su(2) such that = Cay(A^) will solve 
(5.12). The advantage is that su(2) is a linear space, and that no constraints need 
to be imposed on A^ n \ as the range of the Cayley map is automatically contained 
within SU(2). A standard nonlinear solver can therefore be used to find A^ n \ 
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Computational savings. Significant computational savings can be obtained by 
rewriting the Cayley map in a more convenient form. We recall from (2.6) that su(2) 
is isomorphic with M 3 , and we denote the vector representation of by a( n ) £ M 3 . 
A small calculation then shows that the Cayley transform can be expressed as 

[/(«) = Cay(^ (n) ) = 1 ? ( (1 - ||a( n )|| 2 )/ + 2A^ n A , (5.13) 

1 + ||a( n )|| v 1 



so that 



£/■(») - I 



2x2 



^(«) _ ||a (n) || 2 ^ 



1 + ||a( n )|| 2 

The terms proportional to V in (5.12) can then be written as 



2x2 



Re 



2r 


i + 


ll a(n) H 2 




2r 


i + 


ll a(n) H 2 




-2T 


i + 


llaWII 2 



Re[(^a a (A^ - ||aW|| 2 J 2x2 )^ 

^— -2 (Re [(^jV^'Vl " l|a (n) || 2 Re f^ ft )V„^l) 

1 + ||a( n )|| V L J L J/ 

(a^ xx( n ) + ||a( n )|| 2 x( n )) Q , (5.14) 



where we have used the expression (2.11) for the Hopf fibration to write = 
(ip n )^a a tp n , as well as the identity 

3 
/3=1 

= i(aW) a -(aWxxW) ai 

which follows easily from (2.7). 

Similarly, the terms in (5.12) involving the derivatives of the Hamiltonian can be 
written using (4.13) as 



Re(V)t(i<7 a )n t ff(^+V2) 



1 



x n x V< 1/2 - (a™ • x n )VF™ +1/2 - (a™ x x") x VH™ +1/2 ) , 



1 + ||a n | 

where H S 2 is the original point vortex Hamiltonian (3.2). 

Combining this expression with (5.14), we get that the first-order equations (5.12) 
for £/"(") are equivalent to the following nonlinear equation for a^ n ': 

4r 



h 



(a W x X W + ||aW|| 2 x( n )) 



X W x Vtf™ +1/2 - (aW • x("))Vtf™ +1/2 - ( a W x X W) x Vtf™ +1/2 . (5.15) 
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While this equation looks considerably more complicated than the original equation 
(5.12), it can in fact be reduced to a system of easier equations, described in the 
following proposition. To set the notation, we decompose sS n ^ into a part a^ 
perpendicular to x( n ), and a part //x( n ) along x( n ): 

a (n) =a W +/ix (n)_ ( 51 g) 

Proposition 5.2. There are two distinct sets of solutions for the direct first-order 
equation (5.15), depending on the value of fx: 

1. For the first class of solutions, n = and a^ n ) = a_j^ satisfies 
AY 

= x ( X W x VH n+1 / 2 ) + (xW • V# n+1 / 2 )x(") x aK (5.17) 



2. For the second class of solutions, 
h 

if 



ix = --(xW • VH n+l ' 2 + af } • (VH n+l / 2 x X W)). (5.18) 



Substituting (5.16) and (5.18) into (5.15) yields an implicit equation for aj™ . 

Proof. To simplify the notation, we omit the superscripts '(n)' from the quantities 
involved, and we denote the gradient of the Hamiltonian simply by VH. In this 
way, (5.15) becomes 

4r 

— (a x x + ||a|| 2 x) = x x VH - (a • x)V# - (a x x) x VH, (5.19) 

to be solved for a. We substitute the decomposition (5.16) into this equation to 
obtain 

4r 

— (a ± x x + ||a|| 2 x) = x x VH - fiVH - (a± x x) x VH, (5.20) 
where we have used the fact that x 2 = 1. Taking the dot product with x gives 

4r 9 

_ || a || 2 = -/i(x- VH) -& ± ■ VH. (5.21) 
Taking the cross product of (5.20) with x results in 

4r 

- — a ± = (x x VH) x x - fi(VH x x) + (x • VH)(a ± x x). (5.22) 

h 

We now take the dot product of this equation with aj_ to get 
4r o 

— r- ll a -Lll = a ± • ^ H - ^ a -L • (ViJ x x), 
h 

and we subtract this equation from (5.21) to obtain the following quadratic equation 
for fj,: 

4r 

—I? + /x(x • VH + a ± • (VH x x)) = 0, 

the solutions of which are [i = and /U given by (5.18). When fi = 0, a = a^, and 
(5.22) simplifies to (5.17). ■ 
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We stress that the solution for /i / is as valid from a physical point of view 
as the solution for [i = 0, and there is a priori no reason to choose one over the 
other. In fact, we will see that consistently choosing the solution with results 
in a different numerical algorithm, which we have termed the fi-Hopf integrator, 
which has the same good properties (symplecticity, second-order accuracy) as the 
integrator obtained by choosing /x = 0. We will present a numerical comparison 
between both integrators at the end of Section 6.1 but for now we choose the solution 
corresponding to /i = in proposition 5.2, as it gives rise to a simpler equation for 
a (»). 



The adjoint equations. Having reduced the direct first-order equations (5.10) 
to a nonlinear equation on the Lie algebra 0u(2) , we now turn to the adjoint 
first-order equations (5.11). For these equations, we write the update map as 

<pM = [/("^V"" 1 ), (5.23) 



so that 



p(n-l) = (y(n-l))t ¥ ,(n) j 



and we substitute this into (5.11) to get 



Re 



0, (5.24) 



which is very similar to (5.12): to obtain (5.24) from (5.12), we just have to make 
the substitution 

[/(«) ^ ( [ /(«-i))t ) h^-h. 

Using the Cayley map, we now write 

= Cay^- 1 )), so that (U^Y = CayC-y^"" 1 )), 

where the matrix A^ n ~ l \ or rather its vector form a( n_1 ), satisfies the following 
equation: 

_ ^(-a^-D x X W + lla^fxW) = 
h 

xW x Vtf § T 1/2 + (a (n " 1} • xW)Vfl£- 1/2 + (a^" 1 ) x X W) x VH^ l/ \ (5.25) 

which was obtained from (5.15) by making the substitution a^ n ) i— > — a^ n_1 ^ and 
h i — y — h. Note that in the above, x( n ) depends on a^ n_1 ^: 

X W = 7r(</>)) = vr = tt (CayiA^- 1 ^^- 1 ^ , (5.26) 

where y>( n_1 ) is given and A^^ 1 ^ G su(2) is the matrix form of a( n_1 ) S M 3 . The 
adjoint equations will therefore be computationally somewhat more expensive to 
solve than the direct equations, since now x^ n ) also depends on the Lie algebra 
element a^" -1 ) . 

The solutions of the equations (5.25) can now be described by the following propo- 
sition: 
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Proposition 5.3. There are two distinct sets of solutions for the adjoint first- order 
equation (5.25), depending on the value of fi: 

1. For the first class of solutions, fi = and &( n ~^ = aj™ ^ satisfies 
AV 

__ a (n-l) = x (n) x ( x (n) x V ffn-l/2^ _ ^(n) . VH n-l/2^(n) x 

For £/ie second class of solutions, 

/x = A( X (») . vff"- 1 / 2 - af _1) • (V^™- 1 / 2 x xW)). 

Substituting this and a( n_1 ) = a^™ ^ + /xx( n ) into (5.25) yields an implicit 

,. f (n-l) 
equation for aj_ 

In 6ot/i cases, x( n ) depends on a( n_1 ) and is found from x^™ -1 ) using (5.26). 

Again, we will focus from this point on the solution with /i = 0, even though the 
other solution is also physically admissible. See the last part of Section 6.1 for a 
more in-depth discussion. 



5.4 Computational algorithm 

We now summarize the developments of the previous sections to arrive at the follow- 
ing numerical algorithm. We restore the index k = 1, . . . , N labeling the individual 
vortices. Given initial conditions V 9 ^ -1 ; & = 1, • • • for the point vortices, apply 
the following algorithm: 



1. For each point vortex, find a^ n ^ G M 3 by solving the adjoint equation: 



4T 
— —a 



(n-l) 



/( -xl n) x(xl n) xV^ +1 / 2 )-(xi" ) -VF n+1 / 2 )xi n) xai n - 1) . (5.27) 

Here x[ depends on aCu ^ and X ' through (5.26). 

2. Let 4 n) = Cay(al n - 1 Vl n_1) and put x< n > = vr^). 

3. For each point vortex, find al E M 3 by solving the direct equation: 

^af = x[ B) x (xf x Vr +1 /2) + ( X W . V ^ +1 / 2 )x( n) x a< n) . (5.28) 

4. Let = Cay(al n Vl n) and put xi™ +1) = 

5. Repeat. 
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As we shall see below, the resulting flow is second-order accurate, symplectic, and 
preserves the unit-length constraint by construction. The advantage of using this 
version of the discrete equations of motion rather than the original formulations (5.4) 
or (5.7) is that the equations (5.27) and (5.28) for a( n_1 ) and a( n ) are defined on 
su(2) N , a linear space, and hence can be solved using a standard nonlinear solver. We 
stress, however, that these equations are equivalent to the original discrete equations 
of motion defined directly on (S 3 )^. 

Properties of the variational integrator. Since we have started from the mid- 
point discretization of a continuous Hamiltonian, the resulting integrator will be 
second-order accurate, symplectic, and (by construction) unit-length preserving (see 
Marsden and West [2001]). The symplectic form preserved by the numerical algo- 
rithm is not the weighted area form on (S 2 )^ given in (3.1) but the two-form 

where is the discrete Lagrangian given in (5.2), which is an 0(h) perturbation 
of it (see Rowley and Marsden [2002]). 



6 Numerical results 

Throughout this section, we will compare the behavior of the Hopf integrator with 
a number of other integrators: 

1. A standard, non- variational integrator: we chose a standard explicit 4th-order 
Runge-Kutta method (RK4), composed with projection onto the unit sphere 
in order to preserve the unit-length constraint. It is well known (see e.g. 
Haircr, Lubich, and Wanner [2002]) that the resulting method is still 4th- 
order in time. 

2. The implicit midpoint method on § 2 , given by 

n+l v n ! x "+l/2 x x «+l/2 

h 47T ^ 3 1 + n 2 "+1/2 . n+1/2 " ^ ' > 

j^k k X j 

This is just the standard midpoint method, applied to the vortex equations 
(1.6). Note that for this vector field, the implicit midpoint method in fact stays 
on the unit sphere without the need for an explicit projection. To see this, 
take the dot product of both sides of the equation with x^ +1//2 and observe 
that the right-hand side vanishes, so that we get 

l x fc _x fcJ x fc - u ' 
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which is equivalent to ||x^|| = ||x]J + , i.e., the length is preserved. This is 
not a general feature of the midpoint method, but is a consequence of the 
particular form of the point vortex equations. Furthermore, the midpoint 
method is symmetric under the interchange x^ o x^ +1 , h <-> —h, and as a 
result the method seems to have excellent long-term conservation properties. 
To the best of our knowledge, it is not known if this version of the midpoint 
method on the sphere is symplectic. 

3. The Lie group method of Eng0 and Faltinsen [2002], applied to the point 
vortex equations of Section 3. This second-order method preserves the vortex 
moment explicitly, and is a self-adjoint Lie group method, resulting in bounded 
energy error. However, this method is not symplectic (see Zhong and Marsden 
[1988]). This method is implemented by solving 

yx = Ad: xp( „ ?) (yo), with £ = - (VH(y ) + VH( Vl )) , 

where H is the point vortex Hamiltonian (3.2), yo,Ui are the point vortex 
locations in M 3 , viewed as the dual so(3)* of the Lie algebra of the rotation 
group, and Ad*(y ) = gyog' 1 - 

Our conclusion is that all three geometric methods outperform 4th-order Runge- 
Kutta with projection over moderately long time intervals. The Hopf and Lie- 
Poisson algorithms are roughly comparable in performance, while the midpoint 
method generally does better than either one: for one thing, the midpoint method 
conserves the vortex moment, a linear constant of the motion, exactly. 

However, the conservation properties of the midpoint algorithm seem to be some- 
what coincidental, and rely on the fact that for the point vortex equations, the 
algorithm stays on the unit sphere without reprojecting. It is therefore not clear 
how to generalize this algorithm to obtain, for instance, higher-order integrators 
with similar conservation properties. By contrast, for the variational Hopf integra- 
tor it suffices to start from a higher-order version of the discrete Lagrangian (5.1) 
to obtain a higher-order variational Hopf integrator. 

6.1 Stable relative equilibria of vortex rings 

Polvani and Dritschel [1993] have investigated the behavior of a ring of TV equidistant 
vortices with the same strength T, placed on a circle of fixed longitude on the sphere 
(see Figure 6.1). They found that this configuration is a stable relative equilibrium, 
provided that N < 7 and that the longitude is higher than a certain critical value 
(dependent on N). For the case N = 6, the critical longitude is given by z c = 1j\fh 
and the stable relative equilibria satisfy zq > z c . The vortex ring rotates around the 
z-axis with angular velocity £1 = (N — 1) j~ ■ 

For our simulation, we choose N = 6,T = 1/6 and zq = 0.92, so that ~ 0.397 and 
the period T ~ 15.819. The motion of the first vortex over a number of periods is 
illustrated in Figure 6.1. 
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t 

Figure 6.1: Left: Initial conditions for the 6-vortex Polvani-Dritschel vortex 
ring. Right: x, y and z-component of the first vortex in the Polvani-Dritschel 
simulation, where the time-step h = 0.1. The trajectory is clearly seen to be 
periodic. 
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Figure 6.2: Comparison of the energy and momentum preservation between 
the geometric integrators, including our Hopf method, and 4th-order Runge- 
Kutta over moderate integration times. Runge-Kutta exhibits a small secular 
drift in the energy and momentum (left panes), while the energy and momen- 
tum are practically conserved for the geometric integrators (right panes) . Here 
h = 0.1 and a = 0.0. The small drift in the energy-momentum error for the 
Hopf integrator is due to round-off error. 



Comparison with other integrators. We next turn to the energy and momen- 
tum conservation properties of the numerical integrator. We simulate the motion of 
the Polvani-Dritschel vortex ring with time step h = 0.1 and regularization param- 
eter a = 0.0, for T = 1000 units of time, using all four integrators. 

In Figure 6.2 we have plotted the energy error AE := E(t n ) — E(to) (top row of 
figures) and the moment error AM := ||M(i n ) — M(to)|| (bottom row of figures) as 
a function of time. Runge-Kutta exhibits a small drift in the energy and momentum 
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error (left panes), whereas the energy and momentum are practically conserved for 
the geometric integrators (center panes). For the Hopf integrator, there is a much 
smaller drift in the energy and moment error, which is presumably due to the fact 
that one step of the Hopf integrator consists of a step with the forward method (5.28) 
followed by a step with the adjoint method (5.27). The error in these underlying 
one-step methods almost cancels out, up to some round-off error which is seen on 
the figure. 

On the right panes of Figure 6.2, it can be seen that the Hopf integrator exhibits 
a small drift in the energy and moment, which is of the order 10 -13 — 10 -12 over 
10 3 units of time. We conjecture that this is due to the way in which the Hopf 
integrator is implemented. Since a relative equilibrium is being tracked, we compute 
an essentially constant Lie group element at each update, and we multiply the 
previous vortex location c/^™ -1 ) by this matrix to obtain the current location ip( n \ 
As U( n ' is essentially constant, round-off errors in the definition of <p accumulate 
rather than cancelling out on average, resulting in the slight linear trend observed 
in the figures. This roundoff error accumulation could be mitigated by the use of 
compensated summation techniques. 



Numerical order calculation. We know from theoretical considerations that 
the Hopf integrator is second-order accurate, and so are the two other geometric 
methods. We now illustrate this statement by comparing the solution trajectories 
generated by the Hopf integrator with the exact trajectories. For 10 choices of time 
step h between 10 -4 and 10 _1 we run the simulation over T = 100 units of time and 
we compute the absolute error between the numerical and the exact solution. We 
consider only the first vortex, since the trajectories of the other vortices differ from 
the first by a rigid rotation. More precisely, for each integrator we do the following: 
if x exact (i n ) is the exact position of the first vortex at time t n = nh and x™^ 1 is the 
numerical trajectory, then we compute 



Ah := max 



r. \ n,h 
^cxact \tn ) X 



int 



for each of the selected time steps. 

Figure 6.3 (left) shows a plot of absolute errors versus time steps for the three 
geometric integrators. All three integrators are of roughly second order. On the 
right pane of Figure 6.3, we have plotted the obtained accuracy for each of the 
three geometric methods as a function of the expended CPU time. We see that the 
Hopf and Lie-Poisson methods are roughly equivalent in terms of efficiency, with the 
midpoint method slightly more efficient, but we caution against reading too much 
into this figure, as the results depend heavily on the implementation of the method. 



The /i-Hopf integrator. For our last simulation involving the Polvani-Dritschel 
vortex ring, we revisit the choice of the /U-component in propositions 5.2 and 5.3. 
We have so far chosen the solution corresponding to fi = 0, as this solution was 
computationally easier to characterize. There is, however, no a priori reason to 
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Figure 6.3: Left: Absolute error for each of the three integrators for the 
Polvani-Dritschel vortex ring over T — 100 units of time. All three integrators 
are second-order accurate in time. Right: Absolute error as a function of CPU 
time expended, again for T — 100 units of time. All three integrators show the 
same scaling in terms of attaining a certain accuracy for a fixed computational 
time. 




Figure 6.4: Comparison of the regular Hopf and the /i-Hopf algorithm. Left: 
Simulation of the x-component of the first vortex, step size h = 0.5. The solid 
line represents the exact trajectory, while the dashed and the dash-dotted line 
represent the Hopf and the /i-Hopf algorithm, respectively. Both numerical 
algorithms stay close to the exact trajectory over long integration times. Right: 
absolute error versus step size for the two Hopf algorithms. Both algorithms 
exhibit similar second-order accuracy in time. 



prefer the /i = solution over the one for which [i ^ 0. To corroborate this 
statement numerically, we have implemented the Hopf numerical algorithm choosing 
the solution corresponding to fj, ^ for the direct equation (5.15) and the adjoint 
equation (5.25). The resulting numerical algorithm is not equivalent to the original 
Hopf algorithm, and has been termed the fj>-Hopf algorithm. 

In Figure 6.4, we compare both versions of the Hopf algorithm. For the simulation 
on the left, we simulated the motion of the Polvani-Dritschel ring over T = 2000 
units of time with step size h = 0.5. The last 100 units of time of the simulation 
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are shown for both integrators, together with a plot of the exact solution. We see 
that both integrators, although not equivalent, track the exact solution well. On 
the right, we have plotted the absolute error versus the step size employed, and we 
see that both the regular Hopf as well as the /i-Hopf integrator are second-order 
accurate. 

6.2 The spherical von Karman vortex street 

An important class of relative equilibria consists of the single and double von 
Karman vortex streets on the sphere described by Chamoun, Kanso, and New- 
ton [2009]. The single vortex street consists of two staggered arrays of vortices, each 
consisting of ./V equidistant vortices of strength T, at fixed colatitudes <p = (j>\ and 
<f> = it — <f>x , together with vortices of strength T n and T s at the north and the south 
pole, respectively (see Figure 6.5). 




Figure 6.5: The spherical von Karman vortex street consists of two staggered 
latitudinal layers of 5 vortices, together with vortices at the North and South 
pole of the sphere. 

For the simulations in this section, we take the number of vortices in each ring to be 
N = 5, and we set the colatitude equal to <pi = ir/3. The vortex strength for the ring 
vortices is set equal to unity, V = 1, while the polar vortices satisfy T n = —T s = 1/2. 
This configuration forms a relative equilibrium which rotates around the z-axis with 
period T = 10.85. Based on the behavior of the planar Von Karman vortex street, it 
is believed that this relative equilibrium is unstable, although no rigorous stability 
analysis exists, to the best of our knowledge. 

Short integration times. In Figure 6.6 we compare the energy behavior of the 
Hopf variational integrator with that of the 4th-order Runge-Kutta method, for 
short amounts of time. Both simulations use a step size h = 0.1 and no regulariza- 
tion, i.e. a = 0. Over the first 100 units of time, the vortices move as a relative 
equilibrium, and we see that the variational integrator preserves the energy almost 
exactly, in contrast to Runge-Kutta. When the simulation time is close to t = 100, 
a small uptick in the energy and moment error can be seen. This deviation from en- 
ergy/moment preservation signals the onset of instability of the relative equilibrium, 
as will become clear from longer simulations. 
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Figure 6.6: For a relative equilibrium, the variational integrator (line) pre- 
serves the energy (left) and moment (right) almost exactly, compared to 4th- 
order Runge-Kutta (dashed), until the onset of instability close to t — 100. 



Longer integration times. For longer integration times, the equilibrium be- 
comes unstable and breaks up, leading to aperiodic motion of the vortices. In this 
regime, the energy is no longer exactly preserved by the Hopf integrator, but ex- 
hibits bounded oscillations, as is to be expected from a symplectic integrator. In 
Figure 6.7 we have plotted both the energy error as well the error in the norm of the 
vortex moment for each of the integrators. The simulation uses time step h = 0.3 
and regularization parameter a = 0.25 and runs over 1000 time units. 




Figure 6.7: Comparison of the energy (left) and momentum conservation 
(right) between the geometric integrators and 4th-order Runge-Kutta. The 
step size is h — 0.3 and the regularization parameter is a — 0.25. Hopf integra- 
tor: dark blue, solid line; Lie-Poisson integrator: green, dashed line; projected 
midpoint integrator: red, dash-dot line; RK4: light blue, dotted line. 



The geometric integrators clearly exhibit bounded oscillations in the energy and the 
norm of the moment, whereas these quantities grow linearly in time for 4th-order 
Runge-Kutta. 
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6.3 Collapse of three vortices 

It is well known that certain configurations of point vortices on the sphere will 
collapse to a point in finite time. For 3 vortices, necessary and sufficient conditions 
for collapse were given by Kidambi and Newton [1998] while Sakajo [2008] identified 
an open set of initial conditions for collapse of 4 vortices. We focus here on the case 
of 3 vortices. 

We simulate the motion of three vortices with strengths T± = T2 = 1, T3 = —1/2 
placed at the vertices of a triangle with side lengths Z12 = y/S/2, Z23 = v2/2 and 
/31 = 1. For this configuration, it can be calculated that collapse occurs after 
r_ = 47r(v23 — y/l7) = 8.4537 units of time. The trajectories of the colliding 
vortices are shown in Figure 6.8. Note that these initial conditions are for the 
unregularized system, i.e. (1.6) with a = 0. Adding some regularization to the 
system effectively amounts to imposing a minimum distance on the vortices and 
will prevent the vortex configuration from collapsing to a single point. 




Figure 6.8: Trajectories of 3 colliding vortices, for the initial conditions de- 
scribed in the text. 

We simulate the system first with a moderate time step, h = 0.1, and some amount 
of regularization, a = 0.10, for T = 15 units of time. On Figure 6.9 we see that 
the geometric integrators, including the Hopf integrator, perform a better job of 
preserving the energy and vortex moment than Runge-Kutta: while all integrators 
show some buildup in the energy and moment error around the time of the collapse, 
the energy and moment return to their original values after the collapse for the 
geometric integrators, but settle down at a different value for Runge-Kutta. 

For long-term simulations, this effect is even more pronounced. After the near- 
collapse, the three vortices travel past each other and (nearly) collapse again at 
a later time, a situation which repeats itself periodically afterwards. Figure 6.10 
shows that for every collapse event, the Runge-Kutta simulation incurs a jump in 
the energy and the moment, whereas the geometric integrators manage to preserve 
these invariants without any noticeable secular trend. Note also that the period of 
time between two collapse events increases gradually for Runge-Kutta, but stays 
(roughly) constant for the geometric integrators. 
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Figure 6.9: Energy (left) and moment (right) conservation for the geometric 
integrators and 4th-order Runge-Kutta close to vortex collapse, which happens 
for the unregularized system at t ps 8.45. While RK4 conserves the energy and 
moment better than the geometric integrator up to the collapse, the energy 
and moment settle down in a different value after collapse. Here, h = 0.1 and 
a = 0.1. 
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Figure 6.10: Numerical simulation of colliding point vortices for T = 500 
units of time, where h = 0.1 and a = 0.10. At regular instances of time, there 
are collapse events, but whereas the energy for the geometric integrators quickly 
returns to its original value after a collapse event, the energy and moment for 
Runge-Kutta increases in jumps at each collapse. 



6.4 Large ensembles of vortices 

We simulate the motion of 40 vortices of strength T = 1/8, distributed randomly 
over the surface of the sphere. The regularization parameter is set to a = 0.1 and 
the step size is h = 0.01. For short integration times (T = 20), as can be seen from 
Figure 6.11, both the Hopf and the Lie-Poisson integrator exhibit bounded energy 
error, and they preserve the unit length constraint almost up to machine precision. 

For moderate integration times (T = 200, with again h = 0.01 and a = 0.1), we have 
that the energy is almost preserved by all three integrators, up to a small, bounded 
error; see Figure 6.12. For the vortex moment, we have that the midpoint integrator 
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Figure 6.11: Left: Energy error for the motion of 40 random point vortices. 
Right: Preservation of the unit length constraint. 



seems to preserve this invariant exactly, while both the Hopf and the Lie-Poisson 
integrator exhibit a small error consistent with the order of the discretization. In 
this figure, the midpoint method is not included, as it preserves the vortex moment 
exactly. 
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Figure 6.12: Energy error (left) and moment error (right) for the motion of 
40 random point vortices over 200 units of time. The geometric integrators 
exhibit similar conservation properties for the energy and the moment. 



7 Conclusions and outlook 



In this paper, we have used the Hopf fibration to construct a linear Lagrangian 
on the three-sphere S 3 , whose Euler-Lagrange equations project down to the point 
vortex equations on S 2 . We have used this Lagrangian formulation to construct 
a variational integrator for point vortices on the sphere, and we have used the 
isomorphism between S 3 and SU (2) to ensure that the unit-length constraint in § 3 
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is automatically satisfied. Below, we discuss some possibilities for future research. 

Extension to higher-order integrators. Our Lagrangian approach to the con- 
struction of discrete point vortex integrators can be extended without major diffi- 
culties to the construction of integrators whose numerical order is higher than two. 
It suffices to take a discretization in (5.1) which is of higher than second order. 
The standard theory of Lagrangian variational integrators (see Marsden and West 
[2001]) then ensures that the resulting discrete equations of motion will have the 
same order of accuracy as the discrete Lagrangian. 

Other Lie group methods. We have constructed a Lie group variational in- 
tegrator by discretizing the linear Lagrangian (4.5) and writing the update map 
in terms of 5 l ?7(2)-matrices. Another approach is due to Bou-Rabee and Marsden 
[2009] (see also Kobilarov and Marsden [2011]). Here the Lagrangian is first written 
in a left-trivialization of TSU(2) = SU(2) x su(2) by mapping (g, g) to (g, £), where 
g _1 <7 = £, and this equation is then discretized and added to the variational principle 
using a Lagrange multiplier, giving rise to the so-called Hamilton-Pontryagin vari- 
ational principle. A similar variational principle, known as the Clebsch variational 
principle, was pioneered in Cotter and Holm [2009]. It would be of considerable 
interest to discretize our linear Lagrangian (4.5) using these augmented variational 
principles and to compare the properties of the resulting discrete mechanical system 
with our straightforward discretization. 

Statistical mechanics of large numbers of point vortices. The statistical 
theory of vortex motion set forth in Onsager [1949] predicts that, under certain 
energetic conditions, like-signed vortices will tend to cluster over time. We have 
not attempted to extract any statistical information from the simulation of large 
number of vortices on the sphere, but it would be of considerable interest to do so. 
To alleviate the 0(iV 2 )-cost of computing the point vortex Hamiltonian while main- 
taining the symplectic nature of the integrator, a geometric fast multipole method 
like the one developed in Chartier, Darrigrand, and Faou [2010] could be used. 

Vorticity distributions on the sphere and other surfaces. Point vortices 
represent the simplest non-trivial distributions of vortices on the sphere. The meth- 
ods proposed in this paper are expect to generalize without any significant difficulty 
to the case of vortex blobs or patches of vorticity on the sphere (see Chorin [1973]; 
Newton [2001]). 

To treat vortical distributions on other surfaces, the following construction from 
prequantization could be used. Recall that the Hopf fibration is the fiber bundle 
associated to the quantum line bundle on S 2 associated with the area form; see e.g. 
Woodhouse [1992]. For the motion of point vortices on a surface S with integral 
area form, one can follow a similar route and lift the motion of the vortices to (the 
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principal fiber bundle associated to) the quantum line bundle, for which a similar 
relation as (2.14) will continue to hold. 

Moreover, PDEs such as the KdV and nonlinear Schrodinger equation can also be 
formulated using a linear Lagrangian, and we hope that the methods introduced in 
this paper may be useful for these systems as well. 

Acknowledgments 

We would like to dedicate this paper to the memory of Hassan Aref, whose kind 
encouragement and insightful remarks at the 2010 SIAM-SEAS meeting provided 
the initial stimulus for this work. Furthermore, we would like to thank J. D. Brown, 
C. Burnett, F. Gay-Balmaz, M. Gotay, D. Holm, E. Kanso, S. D. Kelly, P. Newton, 
T. Ohsawa, B. Shashikanth and A. Stern for stimulating discussions and helpful 
remarks. 

M. L. and J. V. are partially supported by NSF grants DMS-1010687 and DMS- 
1065972. J. V. is on leave from a Postdoctoral Fellowship of the Research Foundation- 
Flanders (FWO-Vlaanderen). This work is supported by the irses project GE- 
OMECH (nr. 246981) within the 7th European Community Framework Programme. 

A Analysis of a planar vortex integrator 

In this appendix, we show that the integrator of Rowley and Marsden [2002] for 
point vortices in the plane shares a number of remarkable features with the Hopf 
integrator, which stem from the fact that both systems are derivable from a linear 
Lagrangian. Similar observations, but for the numerical integration of canonical 
Hamiltonian systems, were made by Brown [2006]. 

Decomposition into one-step methods. Rowley and Marsden [2002] start from 
the linear Lagrangian (1.3), which they discretize by setting 



where a £ [0, 1] is a real interpolation parameter. The equations of motion derived 
from this Lagrangian are given by 



where z n+a := (1 — a)z n + az n+ \ and f(z) is the right-hand side of the vortex 
equations (1.1). It turns out that for a = 1/2, they can be written as the composition 
of a one-step method and its adjoint. To see this, we specialize to the case a = 1/2 
and use the fact that the original Lagrangian L is linear in the velocities to write 




Zn+2 %n 



af(z n+a ) + (1 - a)f(z n+1 ^ a ), 



(A.l) 



2h 



Ld{z , Zi) = L(zy 2 , *i ) - L(zi/2, Zq) 



A Analysis of a planar vortex integrator 



34 



and we define L d)+ (z , zi,h) := L(z 1 / 2 ,z 1 ) and L d _(zo> *i, h) := -L(z 1/2 ,z ), so 
that L d = L d + + L d _. Consider the adjoint L* d of a discrete Lagrangian L d , which 
is defined by L* d (zQ, zi, h) := —L d (zi, zq, —h) (see Marsden and West [2001]). Then, 
we have that 

L* d)+ (z ,zi,h) = L d ,-(zo,zi,h), 

and vice versa. This definition is motivated by the fact that the adjoint of the 
discrete Euler-Lagrange flow of a discrete Lagrangian is given by the discrete Euler- 
Lagrange flow of the adjoint discrete Lagrangian. 

The composition of the discrete Euler-Lagrange flow of two discrete Lagrangians is 
given by the discrete Euler-Lagrange flow of a composition discrete Lagrangian that 
is the sum of the two original discrete Lagrangians. As a result, the discrete Euler- 
Lagrange flow for L d is given by the composition of the discrete Euler-Lagrange 
flows for L d<+ and its adjoint L* d , = L d - . These discrete flows can be viewed as 
one-step methods, and are typically only first-order accurate, but their composition 
is symmetric and therefore has even order of accuracy, and is, in particular, second- 
order accurate. 

Lastly, we remark that for the point-vortex Lagrangian (1.3) the discrete Lagrangians 
L dt+ and L d ~ coincide, which means that each of them are self-adjoint. As a result, 
the underlying one-step method is second-order. In fact, it can easily be seen that 
for a = 1/2, the point vortex equations (A.l) can be written as the composition of 
the implicit midpoint method 

■2-n+l Z n . . 

T — J\ z n+l/2) 

with itself. This method is clearly second-order accurate. 

For the case of point vortices on the sphere the Lagrangians L d + and L d _ still 
coincide, but in order to recover the equations of motion (5.4) and to enforce the 
constraint (</? n+1 , y? n+1 ) = 1, different constraint forces have to added to the discrete 
flow. As a result, the underlying one-step methods (5.10) and (5.11) no longer 
coincide and are individually only first-order accurate, although their composition 
is second-order accurate. 



The choice a = 0, 1 for the interpolation parameter. The method (A.l) is 
implicit for all choices of a except a = 0, 1, in which case the equations become 



2h 



f(z, 



n+lj 



This method turns out to be the symmetric explicit midpoint method (see Hairer, 
Lubich, and Wanner [2002]), which is well-known to exhibit parasitic oscillatory 
solutions. These solutions can easily be observed in the dynamics of point vor- 
tices: in Figure A.l, we have plotted the energy error for a simulation of a four- 
vortex problem with vortex strengths V = (.1, .3, — .2, — .4) and initial conditions 
z = (0, .5i, 1, .7+ M). For the simulation where a = 0.9 the energy error is bounded, 
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while for the simulation employing a = 1.0 there is a clear linear drift in the energy 
error. The time step used for both simulations was h = 0.1. 

This is in clear contrast to the construction of variational integrators for nondegen- 
erate Lagrangians, for which any choice of interpolation parameter a will result in 
a stable, second-order integrator. 

Similar instabilities exist for the case of point vortices on the sphere: the discrete 
equations (5.6), for instance, exhibit a clear secular drift in the energy, despite being 
variational. 




Time 

Figure A.l: For the four-vortex problem described in the text, the energy 
error exhibits a linear drift for the integrator with a = 1 (solid line) but stays 
bounded whenever a ^ 1; here a = 0.9 is shown (dashed line). 
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